Gardens of forking data
Workspace setup:
Recall from last lecture, the probability that a specific hypothesis, H, is true given a set of data, D, is defined as:
\[ P(H|D) = \frac{P(H)P(D|H)}{P(D)} \]
You can get there via some calculus, but there are other, more intuitive ways of getting to this value. Today, we’re going to build up the intuition of this formula.
Suppose you have a globe representing the planet, and you want to estimate how much of the surface is covered in water. You adopt the following strategy: You toss the globe in the air, and when you catch it, you write down whether the surface under your right thumb is water or land.
Write a function to simulate tosses of this globe. Make sure your function allows you to vary N, the number of tosses, and p, the true proportion of water on the globe.
# function to toss a globe covered p by water N times
sim_globe = function( p=0.7 , N=9 ){
sample(
x = c("W", "L"), # possible values
size = N, # how many draws
prob = c(p, 1-p), # probability of each possibility
replace = TRUE # the same value can be drawn multiple times
)
}
sim_globe()[1] "W" "W" "W" "W" "W" "W" "W" "L" "L"
\[ P(H|D) = \frac{P(H)P(D|H)}{P(D)} \]
What probability distribution would represent the likelihood of a specific sample of water and land given a known propotion of water?
–
The binomial distribution.
Write a function that computes the posterior distribution of \(p\), the true proportion of water, based on a given sample and a flat prior.
\[ P(H|D) = \frac{P(H)P(D|H)}{P(D)} \]
Testing it out:
[1] "W" "W" "L" "W" "W" "L" "W" "W" "L"
[1] 0.000000e+00 8.741894e-12 5.425284e-10 5.990575e-09 3.261806e-08
[6] 1.205399e-07 3.485649e-07 8.509010e-07 1.834811e-06 3.598404e-06
[11] 6.547829e-06 1.121325e-05 1.826302e-05 2.851562e-05 4.294893e-05
[16] 6.270652e-05 8.910077e-05 1.236125e-04 1.678871e-04 2.237272e-04
[21] 2.930816e-04 3.780305e-04 4.807680e-04 6.035804e-04 7.488219e-04
[26] 9.188879e-04 1.116184e-03 1.343097e-03 1.601956e-03 1.895001e-03
[31] 2.224345e-03 2.591936e-03 2.999519e-03 3.448599e-03 3.940405e-03
[36] 4.475850e-03 5.055498e-03 5.679534e-03 6.347729e-03 7.059415e-03
[41] 7.813458e-03 8.608241e-03 9.441644e-03 1.031103e-02 1.121325e-02
[46] 1.214461e-02 1.310092e-02 1.407747e-02 1.506903e-02 1.606992e-02
[51] 1.707401e-02 1.807474e-02 1.906518e-02 2.003808e-02 2.098589e-02
[56] 2.190088e-02 2.277513e-02 2.360067e-02 2.436951e-02 2.507375e-02
[61] 2.570565e-02 2.625773e-02 2.672284e-02 2.709431e-02 2.736600e-02
[66] 2.753241e-02 2.758880e-02 2.753126e-02 2.735684e-02 2.706360e-02
[71] 2.665076e-02 2.611870e-02 2.546910e-02 2.470498e-02 2.383075e-02
[76] 2.285223e-02 2.177672e-02 2.061293e-02 1.937103e-02 1.806258e-02
[81] 1.670044e-02 1.529871e-02 1.387258e-02 1.243815e-02 1.101227e-02
[86] 9.612248e-03 8.255589e-03 6.959637e-03 5.741183e-03 4.616016e-03
[91] 3.598404e-03 2.700509e-03 1.931739e-03 1.298012e-03 8.009482e-04
[96] 4.369673e-04 1.962992e-04 6.189388e-05 8.227801e-06 0.000000e+00
Compute and plot the grid approximate posterior distribution for each of the following sets of observations. In each case, assume a uniform prior for \(p\):
Now assume a prior for \(p\) that is equal to zero when \(p < 0.5\) and is a positive constant when \(p ≥ 0.5\). Again compute and plot the grid approximate posterior distribution for each of the sets of observations in the prior problem.
Suppose the globe tossing data had turned out to be 8 water in 15 tosses. Construct the posterior distribution, using grid approximation. Use the same flat prior as before.
Draw 10,000 samples from the grid approximation from above. Then use the samples to calculate the 90% HPDI for \(p\).
p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior <- rep( 1 , 1000 )
likelihood <- dbinom( 8 , size=15 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
set.seed(100)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )
rethinking::HPDI(samples, prob = .90) |0.9 0.9|
0.3343343 0.7217217
Construct a posterior predictive check for this model and data. This means simulate the distribution of samples, averaging over the posterior uncertainty in \(p\). What is the probability of observing 8 water in 15 tosses?
Now use a prior that is zero below \(p <.5\) and a positive constant otherwise. Repeat each problem above and compare the inferences. What difference does the better prior make?
p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior <- ifelse(p_grid < .5, 0, 1)
likelihood <- dbinom( 8 , size=15 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
set.seed(100)
samples <- sample( p_grid , prob=posterior , size=1e4 , replace=TRUE )
rethinking::HPDI(samples, prob = .90) |0.9 0.9|
0.5005005 0.7097097
dummy_w
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
2 7 35 132 347 684 1081 1630 1776 1667 1273 769 388 184 25
[1] 0.163